Changes in Fv/Fm in elevated nutrients and heat stress
Ana Palacio
2020-12-20
- General project set-up
- Data exploration
- Figure 3: Treatment differences by species and dominant symbiont
- YII GLMs
- Packages used
General project set-up
# Load libraries and sources required to run the script
library(tidyverse)
library(ggthemes)
library(lmerTest)
library(emmeans)
library(multcomp)
library(effects)
library(gridExtra)
library(rstatix)
# Default ggplot settings
Fill.colour<-scale_colour_manual(values = c("black", "gray70", "gray35"))
ggthe_bw<-theme(plot.background=element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
legend.box.background = element_rect(),
panel.background =element_rect(fill = NA, color = "black")
)+
theme_bw()Data exploration
1. Get the files with all the YII by species
YII.data<-read.csv("YII_Data/All_YII_data.csv", header = T)
summary(YII.data)## Sample Date Spp Fragment Treatment
## Ac_288_T21: 2 2017-11-16: 354 Ac:2281 Ac_108 : 24 A :2096
## Ac_101_T0 : 1 2017-11-23: 354 Of:1051 Ac_116 : 24 N :1974
## Ac_101_T1 : 1 2017-11-29: 354 Ss:2687 Ac_119 : 24 N+P:1949
## Ac_101_T10: 1 2017-12-06: 354 Ac_122 : 24
## Ac_101_T11: 1 2017-12-13: 354 Ac_143 : 24
## Ac_101_T12: 1 2018-01-03: 354 Ac_152 : 24
## (Other) :6012 (Other) :3895 (Other):5875
## Replicate YII Genotype Days
## R1:3061 Min. :0.0000 G_62 : 577 Min. :-112.00
## R2:2958 1st Qu.:0.4420 G_48 : 570 1st Qu.: 14.00
## Median :0.5030 G_07 : 462 Median : 65.00
## Mean :0.4988 Ss_20 : 402 Mean : 52.42
## 3rd Qu.:0.5720 Ss_23 : 402 3rd Qu.: 92.00
## Max. :0.6870 Ss_27 : 400 Max. : 156.00
## (Other):3206
## Time_Point Phase TotalSH logSH
## T10 : 354 Baseline : 954 Min. :0.000 Min. :-5.473
## T5 : 354 Heat :1520 1st Qu.:0.020 1st Qu.:-1.709
## T6 : 354 Nutrients:2818 Median :0.057 Median :-1.243
## T7 : 354 Ramping : 514 Mean :0.116 Mean :-1.317
## T8 : 354 Recovery : 213 3rd Qu.:0.158 3rd Qu.:-0.802
## T9 : 354 Max. :1.394 Max. : 0.144
## (Other):3895 NA's :4657 NA's :4657
## D.Prp Community InitialCommunity
## Min. :0.0000 A :2281 A :2281
## 1st Qu.:0.0000 B : 195 B : 186
## Median :0.0000 C1: 538 C1: 468
## Mean :0.1705 C3: 727 C3: 762
## 3rd Qu.:0.0000 D :2278 D :2322
## Max. :1.0000
## NA's :2836Merge/Transform
# Organize data type
YII.data$Date<-as.Date(YII.data$Date, "%Y-%m-%d")
YII.data$Days<-(as.numeric(YII.data$Date) -17485)
#Time as a factor, not as int
str(YII.data)## 'data.frame': 6019 obs. of 16 variables:
## $ Sample : Factor w/ 6018 levels "Ac_101_T0","Ac_101_T1",..: 1 2 10 11 12 13 14 15 16 17 ...
## $ Date : Date, format: "2017-07-26" "2017-08-30" ...
## $ Spp : Factor w/ 3 levels "Ac","Of","Ss": 1 1 1 1 1 1 1 1 1 1 ...
## $ Fragment : Factor w/ 354 levels "Ac_101","Ac_102",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment : Factor w/ 3 levels "A","N","N+P": 3 3 3 3 3 3 3 3 3 3 ...
## $ Replicate : Factor w/ 2 levels "R1","R2": 1 1 1 1 1 1 1 1 1 1 ...
## $ YII : num 0.644 0.576 0.563 0.568 0.645 0.589 0.595 0.606 0.605 0.606 ...
## $ Genotype : Factor w/ 17 levels "G_07","G_08",..: 5 5 5 5 5 5 5 5 5 5 ...
## $ Days : num -112 -77 -36 -27 -9 1 8 14 21 28 ...
## $ Time_Point : Factor w/ 25 levels "T0","T1","T10",..: 1 2 12 19 20 21 22 23 24 25 ...
## $ Phase : Factor w/ 5 levels "Baseline","Heat",..: 1 1 1 1 1 1 3 3 3 3 ...
## $ TotalSH : num NA NA NA NA NA ...
## $ logSH : num NA NA NA NA NA ...
## $ D.Prp : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Community : Factor w/ 5 levels "A","B","C1","C3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ InitialCommunity: Factor w/ 5 levels "A","B","C1","C3",..: 1 1 1 1 1 1 1 1 1 1 ... YII.data$DaysF<-as.factor(YII.data$Days)
YII.data$Spp <- as.factor(YII.data$Spp)
YII.data$Treatment <- as.factor(YII.data$Treatment)
YII.data$Genotype<-factor(as.character(YII.data$Genotype),
levels=c("G_48", "G_62","G_31",
"G_08","G_07", "G_50",
"Of_34","Of_20","Of_6", "Of_31",
"Ss_22","Ss_23","Ss_27", "Ss_28",
"Ss_20", "Ss_24","Ss_30"
)) # D dominance order
YII.data$Community <- factor(YII.data$Community,
levels = c("A", "B", "C3", "C1", "D"))
# Check the data
str(YII.data)## 'data.frame': 6019 obs. of 17 variables:
## $ Sample : Factor w/ 6018 levels "Ac_101_T0","Ac_101_T1",..: 1 2 10 11 12 13 14 15 16 17 ...
## $ Date : Date, format: "2017-07-26" "2017-08-30" ...
## $ Spp : Factor w/ 3 levels "Ac","Of","Ss": 1 1 1 1 1 1 1 1 1 1 ...
## $ Fragment : Factor w/ 354 levels "Ac_101","Ac_102",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment : Factor w/ 3 levels "A","N","N+P": 3 3 3 3 3 3 3 3 3 3 ...
## $ Replicate : Factor w/ 2 levels "R1","R2": 1 1 1 1 1 1 1 1 1 1 ...
## $ YII : num 0.644 0.576 0.563 0.568 0.645 0.589 0.595 0.606 0.605 0.606 ...
## $ Genotype : Factor w/ 17 levels "G_48","G_62",..: 6 6 6 6 6 6 6 6 6 6 ...
## $ Days : num -112 -77 -36 -27 -9 1 8 14 21 28 ...
## $ Time_Point : Factor w/ 25 levels "T0","T1","T10",..: 1 2 12 19 20 21 22 23 24 25 ...
## $ Phase : Factor w/ 5 levels "Baseline","Heat",..: 1 1 1 1 1 1 3 3 3 3 ...
## $ TotalSH : num NA NA NA NA NA ...
## $ logSH : num NA NA NA NA NA ...
## $ D.Prp : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Community : Factor w/ 5 levels "A","B","C3","C1",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ InitialCommunity: Factor w/ 5 levels "A","B","C1","C3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ DaysF : Factor w/ 26 levels "-112","-77","-36",..: 1 2 3 4 5 6 7 8 9 10 ... summary(YII.data)## Sample Date Spp Fragment
## Ac_288_T21: 2 Min. :2017-07-26 Ac:2281 Ac_108 : 24
## Ac_101_T0 : 1 1st Qu.:2017-11-29 Of:1051 Ac_116 : 24
## Ac_101_T1 : 1 Median :2018-01-19 Ss:2687 Ac_119 : 24
## Ac_101_T10: 1 Mean :2018-01-06 Ac_122 : 24
## Ac_101_T11: 1 3rd Qu.:2018-02-15 Ac_143 : 24
## Ac_101_T12: 1 Max. :2018-04-20 Ac_152 : 24
## (Other) :6012 (Other):5875
## Treatment Replicate YII Genotype Days
## A :2096 R1:3061 Min. :0.0000 G_62 : 577 Min. :-112.00
## N :1974 R2:2958 1st Qu.:0.4420 G_48 : 570 1st Qu.: 14.00
## N+P:1949 Median :0.5030 G_07 : 462 Median : 65.00
## Mean :0.4988 Ss_23 : 402 Mean : 52.42
## 3rd Qu.:0.5720 Ss_20 : 402 3rd Qu.: 92.00
## Max. :0.6870 Ss_27 : 400 Max. : 156.00
## (Other):3206
## Time_Point Phase TotalSH logSH
## T10 : 354 Baseline : 954 Min. :0.000 Min. :-5.473
## T5 : 354 Heat :1520 1st Qu.:0.020 1st Qu.:-1.709
## T6 : 354 Nutrients:2818 Median :0.057 Median :-1.243
## T7 : 354 Ramping : 514 Mean :0.116 Mean :-1.317
## T8 : 354 Recovery : 213 3rd Qu.:0.158 3rd Qu.:-0.802
## T9 : 354 Max. :1.394 Max. : 0.144
## (Other):3895 NA's :4657 NA's :4657
## D.Prp Community InitialCommunity DaysF
## Min. :0.0000 A :2281 A :2281 1 : 354
## 1st Qu.:0.0000 B : 195 B : 186 8 : 354
## Median :0.0000 C3: 727 C1: 468 14 : 354
## Mean :0.1705 C1: 538 C3: 762 21 : 354
## 3rd Qu.:0.0000 D :2278 D :2322 28 : 354
## Max. :1.0000 49 : 354
## NA's :2836 (Other):3895Remove / subset timepoints
# Remove baseline values
YII.data<-subset(YII.data, Days>-1)
# Remove recovery values
YII.data<-subset(YII.data, Days<112)
# write.csv(YII.data, "Outputs/Experiment_YII_data.csv", row.names = F)
# YII.Wide<- reshape(YII.data, idvar = "Fragment", timevar = "Days", direction = "wide")
Spp.fragments<-YII.data %>%
group_by(Spp, Genotype, Treatment, Replicate) %>% count(Fragment)
Spp.fragments #write.csv(Spp.fragments, "Outputs/Meassurments_perFragments.csv", row.names = F)
# Subset data
YII.nutrients<-subset(YII.data, Days<80)
YII.heat<-subset(YII.data, Days>75)2. Exploratory graphs
All time points (nutrients + heat stress)
- Colony (Genotype) differences
YII_Colony<- ggplot(YII.data, aes (Days, YII, colour=Genotype)) +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.5)+
stat_summary(fun.y=mean, geom="line", alpha=0.6) + theme_bw()
YII_Colony + ylim(0.0, 0.65) + facet_grid (Spp~Treatment)
YII_Frag_Gen<- ggplot(YII.data, aes (Days, YII, colour=Genotype, group=Fragment)) +
stat_summary(fun.y=mean, geom="line", alpha=0.5) +
theme_bw() + theme(legend.position = "bottom",
legend.title = element_blank())
YII_Frag_Gen + ylim(0.0, 0.65) + facet_grid (Spp~Treatment)
Figure 3: Treatment differences by species and dominant symbiont
DropPlot<-YII.data[which(YII.data$Spp !="Ac" &
YII.data$DaysF==110),]
YII.datab<-YII.data[!(YII.data$Sample) %in% (DropPlot$Sample),]
YII_Treat<- ggplot(YII.datab, aes (Days, YII, colour=Treatment)) +
#stat_summary(fun.data = "mean_cl_boot",geom = "errorbar",
# width = 0.2, position = position_dodge(1) ) +
ggthe_bw+ Fill.colour +
# stat_summary(fun.y=mean, geom="point") +
theme(legend.position="bottom",
strip.background = element_rect(fill="white"))+
scale_y_continuous(limits = c(0, 0.7),
breaks = seq(0, 0.6, 0.1),
expand = c(0, 0),
name=("Fv/Fm")) +
scale_x_continuous(name="Days in the experiment",
limits = c(-1,113),
breaks = seq(0, 113, 15),
expand = c(0, 0))+
annotate("segment", x = 2, xend = 91, y = 0.05, yend = 0.05,
colour = "gray90", linetype=2)+
annotate("segment", x = 79, xend = 91, y = 0.05, yend = 0.20,
colour = "gray90", linetype=3)+
annotate("segment", x = 91, xend = 110, y = 0.20, yend = 0.20,
colour = "gray90", linetype=3)+
annotate("text", x = 45, y = 0.12, label = "Nutrients", size=3)+
annotate("text", x = 99, y = 0.12, label = "Heat", size=3)
#YII_Treat + facet_grid (~Spp) + geom_smooth(span=0.5)
All_SppFif<-YII_Treat + facet_wrap (Spp~Community) + geom_smooth(span=0.5)
All_SppFif
#ggsave(file="Outputs/All_spp_YII_Treatb.svg", plot=All_SppFif, width=7, height=6)Figure 3b (No smooth)
YII_Treat_BW<- ggplot(data=YII.data, aes (Days, YII, colour=factor(Treatment), shape=factor(Treatment))) +
ggthe_bw + Fill.colour+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position = position_dodge(1) )+
stat_summary(fun.y=mean, geom="line", position = position_dodge(1),
linetype=1, alpha=1) +
stat_summary(fun.y=mean, geom="point", size = 2,
position=position_dodge(width=1), alpha=0.8) +
theme(legend.position=c(0.1, 0.8),
legend.title = element_blank(),
strip.background =element_rect(fill=NA)) + # geom_smooth()+
scale_y_continuous(limits = c(0.1, 0.7),
breaks = seq(0.1, 0.6, 0.1),
expand = c(0, 0),
name=expression(~italic("Fv / Fm"))) +
scale_x_continuous(name="Days in the experiment",
limits = c(-1,113),
breaks = seq(0, 113, 15),
expand = c(0, 0))+
annotate("segment", x = 2, xend = 91, y = 0.12, yend = 0.12,
colour = "gray90", linetype=1)+
annotate("segment", x = 79, xend = 91, y = 0.12, yend = 0.20,
colour = "gray90", linetype=1)+
annotate("segment", x = 91, xend = 110, y = 0.20, yend = 0.20,
colour = "gray90", linetype=1)
Figure3b<-YII_Treat_BW + #facet_grid (Spp~.)+
annotate("text", x = 45, y = 0.15, label = "Nutrients", size=3)+
annotate("text", x = 99, y = 0.15, label = "Heat", size=3)+
facet_wrap(Spp~Community)
Figure3b
#ggsave(file="Outputs/Fig_3b_Acer_YII_Treat.svg", plot=Figure3, width=3.5, height=3)Figure 3c (Symbiont based)
YII_Treat_Sym_Spp<- ggplot(data=subset(YII.data), aes (Days, YII, colour=factor(InitialCommunity))) +
ggthe_bw + geom_smooth()+
theme(legend.position="bottom",
legend.title = element_blank(),
strip.background = element_rect(fill="white"))+
scale_y_continuous(limits = c(0.2, 0.7),
breaks = seq(0.1, 0.65, 0.1),
expand = c(0, 0),
name=expression(~italic("Fv / Fm"))) +
scale_x_continuous(name="Days in the experiment",
limits = c(-1,113),
breaks = seq(0, 113, 15),
expand = c(0, 0))+
annotate("segment", x = 2, xend = 91, y = 0.21, yend = 0.21,
colour = "gray90", linetype=2)+
annotate("segment", x = 79, xend = 91, y = 0.21, yend = 0.28,
colour = "gray90", linetype=3)+
annotate("segment", x = 91, xend = 110, y = 0.28, yend = 0.28,
colour = "gray90", linetype=3)
Figure3c<-YII_Treat_Sym_Spp + facet_grid (Spp~Treatment)+
annotate("text", x = 45, y = 0.25, label = "Nutrients", size=3)+
annotate("text", x = 99, y = 0.25, label = "Heat", size=3)
Figure3c
YII GLMs
All spp pooled (does not consider domminant symbiont)
# All spp together
LME1<-lmer(YII ~ Treatment * DaysF * Spp + (1|Genotype),
REML=TRUE, data=YII.data, na.action=na.omit)
lmerTest::step (LME1)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 154 8708.5 -17109
## (1 | Genotype) 0 153 8299.7 -16293 817.54 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## Treatment:DaysF:Spp 0 1.2893 0.020466 63 4851.1 14.11
## Pr(>F)
## Treatment:DaysF:Spp < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment * DaysF * Spp + (1 | Genotype) drop1(LME1, test = "Chisq") #summary(LME1)
anova(LME1) ranova(LME1) # EMMs
All.YII.emm<-emmeans(LME1, ~Treatment * DaysF* Spp)
emmip(LME1, ~DaysF|Treatment|Spp, CIs = TRUE) + theme_bw() + facet_grid(Spp~Treatment)
# Spp responded differently, do separate analysis for each oneA.cer model
Nutrient treatment
Subset Acervicornis data
YII.Acer<-subset(YII.data, Spp=="Ac")
YII.Acer$Nutrients<-"Nutrients"
YII.Acer$Nutrients[YII.Acer$Treatment=="A"]<-"Ambient"Find best model
# 1. Find the best model
YII.Acer$DaysF<-as.factor(YII.Acer$Days)
LME_Acer<-lmerTest::lmer(YII ~ Treatment * DaysF +
(1|Genotype) + (1|Replicate) + (1|Fragment),
data=YII.Acer, na.action=na.omit)
#summary(LME_Acer)
Step.LME_Acer<-step (LME_Acer) # Replicate is not significant
anova(LME_Acer) ranova(LME_Acer)# Replicate is not significant # Drop (1|Replicate)
final_fm <- get_model(Step.LME_Acer)
#summary(final_fm)
LME_Acer1<-lmerTest::lmer(YII ~ Treatment * DaysF +
(1|Genotype/Fragment),
data=subset(YII.data, Spp=="Ac"), na.action=na.omit)
ranova(LME_Acer1) LME_Acer2<-lmerTest::lmer(YII ~ Treatment * DaysF +
(1|Genotype) + (1|Fragment),
data=subset(YII.data, Spp=="Ac"), na.action=na.omit)
ranova(LME_Acer2) anova(LME_Acer1, LME_Acer2) # LME_Acer1 and LME_Acer2 are the same #2. Extract EMMs
Acer.YII.emm<-emmeans(LME_Acer1, ~Treatment | DaysF)
contrast(Acer.YII.emm, "tukey")## DaysF = 1:
## contrast estimate SE df t.ratio p.value
## A - N 0.01123 0.00566 1023 1.985 0.1164
## A - N+P 0.00684 0.00569 1022 1.202 0.4524
## N - N+P -0.00439 0.00562 1023 -0.781 0.7147
##
## DaysF = 8:
## contrast estimate SE df t.ratio p.value
## A - N -0.01626 0.00566 1023 -2.874 0.0115
## A - N+P -0.01497 0.00569 1022 -2.630 0.0235
## N - N+P 0.00129 0.00562 1023 0.229 0.9715
##
## DaysF = 14:
## contrast estimate SE df t.ratio p.value
## A - N 0.00304 0.00566 1023 0.537 0.8533
## A - N+P -0.00902 0.00569 1022 -1.584 0.2528
## N - N+P -0.01206 0.00562 1023 -2.145 0.0815
##
## DaysF = 21:
## contrast estimate SE df t.ratio p.value
## A - N -0.03309 0.00566 1023 -5.849 <.0001
## A - N+P -0.02066 0.00569 1022 -3.629 0.0009
## N - N+P 0.01243 0.00562 1023 2.211 0.0698
##
## DaysF = 28:
## contrast estimate SE df t.ratio p.value
## A - N -0.03049 0.00566 1023 -5.389 <.0001
## A - N+P -0.03612 0.00569 1022 -6.344 <.0001
## N - N+P -0.00562 0.00562 1023 -1.000 0.5768
##
## DaysF = 49:
## contrast estimate SE df t.ratio p.value
## A - N -0.05631 0.00566 1023 -9.952 <.0001
## A - N+P -0.05586 0.00569 1022 -9.811 <.0001
## N - N+P 0.00045 0.00562 1023 0.080 0.9965
##
## DaysF = 65:
## contrast estimate SE df t.ratio p.value
## A - N -0.03088 0.00566 1023 -5.457 <.0001
## A - N+P -0.04943 0.00573 1031 -8.633 <.0001
## N - N+P -0.01855 0.00565 1032 -3.281 0.0031
##
## DaysF = 71:
## contrast estimate SE df t.ratio p.value
## A - N -0.04150 0.00575 1050 -7.212 <.0001
## A - N+P -0.02138 0.00573 1031 -3.734 0.0006
## N - N+P 0.02012 0.00575 1059 3.500 0.0014
##
## DaysF = 76:
## contrast estimate SE df t.ratio p.value
## A - N -0.03999 0.00579 1061 -6.909 <.0001
## A - N+P -0.03217 0.00573 1031 -5.619 <.0001
## N - N+P 0.00782 0.00578 1070 1.353 0.3664
##
## DaysF = 84:
## contrast estimate SE df t.ratio p.value
## A - N -0.02435 0.00644 1226 -3.784 0.0005
## A - N+P -0.03444 0.00634 1194 -5.436 <.0001
## N - N+P -0.01009 0.00652 1252 -1.547 0.2692
##
## DaysF = 89:
## contrast estimate SE df t.ratio p.value
## A - N -0.02714 0.00649 1238 -4.179 0.0001
## A - N+P -0.02383 0.00634 1194 -3.760 0.0005
## N - N+P 0.00331 0.00658 1263 0.503 0.8698
##
## DaysF = 92:
## contrast estimate SE df t.ratio p.value
## A - N 0.01667 0.00676 1295 2.466 0.0367
## A - N+P -0.01762 0.00649 1230 -2.713 0.0186
## N - N+P -0.03429 0.00699 1340 -4.906 <.0001
##
## DaysF = 96:
## contrast estimate SE df t.ratio p.value
## A - N 0.05921 0.00712 1360 8.311 <.0001
## A - N+P 0.02681 0.00712 1352 3.764 0.0005
## N - N+P -0.03240 0.00790 1452 -4.103 0.0001
##
## DaysF = 99:
## contrast estimate SE df t.ratio p.value
## A - N 0.12685 0.00736 1393 17.243 <.0001
## A - N+P 0.04875 0.00736 1388 6.627 <.0001
## N - N+P -0.07810 0.00831 1485 -9.397 <.0001
##
## DaysF = 103:
## contrast estimate SE df t.ratio p.value
## A - N 0.12243 0.00781 1447 15.678 <.0001
## A - N+P 0.10504 0.00749 1406 14.022 <.0001
## N - N+P -0.01739 0.00883 1516 -1.970 0.1200
##
## DaysF = 106:
## contrast estimate SE df t.ratio p.value
## A - N 0.18950 0.00877 1513 21.601 <.0001
## A - N+P 0.14946 0.00821 1478 18.205 <.0001
## N - N+P -0.04004 0.01026 1555 -3.904 0.0003
##
## DaysF = 110:
## contrast estimate SE df t.ratio p.value
## A - N nonEst NA NA NA NA
## A - N+P 0.13704 0.01142 1565 12.001 <.0001
## N - N+P nonEst NA NA NA NA
##
## P value adjustment: tukey method for comparing a family of 3 estimates Acer.YII.emm<-emmeans(LME_Acer1, ~Treatment * DaysF)
# Effect plot options
emmip(LME_Acer, ~DaysF|Treatment, CIs = TRUE) + theme_bw() # interaction plot of predictions
Acer.YII_groups<-cld(Acer.YII.emm, by=NULL) # compact-letter display
Acer.YII_groups<-Acer.YII_groups[order(Acer.YII_groups$Treatment, Acer.YII_groups$Day),]
Acer.YII_groups #write.csv(Acer.YII_groups, "Outputs/Multicomp_AcerYII.csv", row.names = F)Ofav and Ssid
Prepare data sets:
- Time as discrete
#Time as a factor, not as int
str(YII.data)## 'data.frame': 5017 obs. of 17 variables:
## $ Sample : Factor w/ 6018 levels "Ac_101_T0","Ac_101_T1",..: 13 14 15 16 17 3 4 5 6 7 ...
## $ Date : Date, format: "2017-11-16" "2017-11-23" ...
## $ Spp : Factor w/ 3 levels "Ac","Of","Ss": 1 1 1 1 1 1 1 1 1 1 ...
## $ Fragment : Factor w/ 354 levels "Ac_101","Ac_102",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment : Factor w/ 3 levels "A","N","N+P": 3 3 3 3 3 3 3 3 3 3 ...
## $ Replicate : Factor w/ 2 levels "R1","R2": 1 1 1 1 1 1 1 1 1 1 ...
## $ YII : num 0.589 0.595 0.606 0.605 0.606 0.625 0.593 0.567 0.593 0.573 ...
## $ Genotype : Factor w/ 17 levels "G_48","G_62",..: 6 6 6 6 6 6 6 6 6 6 ...
## $ Days : num 1 8 14 21 28 49 65 71 76 84 ...
## $ Time_Point : Factor w/ 25 levels "T0","T1","T10",..: 21 22 23 24 25 3 4 5 6 7 ...
## $ Phase : Factor w/ 5 levels "Baseline","Heat",..: 1 3 3 3 3 3 3 3 3 4 ...
## $ TotalSH : num 0.1262 NA NA NA 0.0401 ...
## $ logSH : num -0.899 NA NA NA -1.397 ...
## $ D.Prp : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Community : Factor w/ 5 levels "A","B","C3","C1",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ InitialCommunity: Factor w/ 5 levels "A","B","C1","C3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ DaysF : Factor w/ 26 levels "-112","-77","-36",..: 6 7 8 9 10 11 12 13 14 15 ... YII.data$DaysF<-as.factor(YII.data$Days)Separate coral species
- Separate phases:
- BaseLine
- C (nutrients only)
- H (Heat challenge)
Ofav:
YII.Ofav<-subset(YII.data, Spp=="Of")
YII.Ofav.0<-subset(YII.Ofav, Days<2)
YII.Ofav.C<-subset(YII.Ofav, Days<77)
YII.Ofav.C.1<-subset(YII.Ofav.C, Days>2)
YII.Ofav.H<-subset(YII.Ofav, Days>75)Ssid
YII.Ssid<-subset(YII.data, Spp=="Ss")
YII.Ssid.0<-subset(YII.Ssid, Days<2)
YII.Ssid.C<-subset(YII.Ssid, Days<77)
YII.Ssid.C.1<-subset(YII.Ssid.C, Days>2)
YII.Ssid.H<-subset(YII.Ssid, Days>75)Baseline: effect of different Symbiodiniaceae taxa before treatments
Ofav (day 0)
LME_Ofav0<-lmer(YII ~ Treatment + InitialCommunity + (1|Replicate),
data=YII.Ofav.0)
#step (LME_Ofav0)
anova(LME_Ofav0) # Treatment not significant ranova(LME_Ofav0) # Replicate is not significant LME_Ofav0<-lm(YII ~ InitialCommunity, data=YII.Ofav.0)
step (LME_Ofav0) ## Start: AIC=-508.95
## YII ~ InitialCommunity
##
## Df Sum of Sq RSS AIC
## <none> 0.057983 -508.95
## - InitialCommunity 1 0.014183 0.072166 -495.19##
## Call:
## lm(formula = YII ~ InitialCommunity, data = YII.Ofav.0)
##
## Coefficients:
## (Intercept) InitialCommunityD
## 0.53931 -0.03376 anova(LME_Ofav0) # EMMs
Ofav.YII.emm0<-emmeans(LME_Ofav0, ~ InitialCommunity)
# contrast(Ssid.YII.emm, "tukey")
# Tukey comparison (do not trust CI to compare EMMs)
plot(emmeans(LME_Ofav0, ~InitialCommunity), comparisons = TRUE) +
coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) +
theme_bw()
Ofav.YII_groups0<-cld(Ofav.YII.emm0, by=NULL) # compact-letter display
Ofav.YII_groups0 # write.csv(Ofav.YII_groups0, "Outputs/Multicomp_OfavYII0.csv", row.names = F)Ssid (day 0)
LME_Ssid0<-lmer(YII ~ Treatment + InitialCommunity + (1|Replicate),
data=YII.Ssid.0)
step (LME_Ssid0) # Treatemnt and replicate are significant :/## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 7 369.23 -724.46
## (1 | Replicate) 0 6 338.52 -665.03 61.431 1 4.586e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## Treatment 0 0.019079 0.009540 2 156 20.809
## InitialCommunity 0 0.087511 0.043756 2 156 95.445
## Pr(>F)
## Treatment 9.754e-09 ***
## InitialCommunity < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment + InitialCommunity + (1 | Replicate) anova(LME_Ssid0) ranova(LME_Ssid0) # EMMs
Ssid.YII.emm0<-emmeans(LME_Ssid0, ~ InitialCommunity)
# contrast(Ssid.YII.emm, "tukey")
# Tukey comparison (do not trust CI to compare EMMs)
plot(emmeans(LME_Ssid0, ~InitialCommunity), comparisons = TRUE) +
coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) +
theme_bw()
Ssid.YII_groups0<-cld(Ssid.YII.emm0, by=NULL) # compact-letter display
Ssid.YII_groups0 # write.csv(Ssid.YII_groups0, "Outputs/Multicomp_SsidYII0.csv", row.names = F)C: Effect of nutrient treatments at control temperature
Ofav
- C Days pooled (days 8-76)
LME_Ofav1.1<-lmer(YII ~ Treatment * InitialCommunity + (1|Fragment),
data=YII.Ofav.C.1)
step (LME_Ofav1.1) # Replicate is significant## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 1046.0 -2075.9
## (1 | Fragment) 0 7 1036.4 -2058.8 19.116 1 1.23e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF
## Treatment:InitialCommunity 1 0.002297 0.001148 2 66
## Treatment 0 0.013198 0.006599 2 68
## InitialCommunity 0 0.173819 0.173819 1 68
## F value Pr(>F)
## Treatment:InitialCommunity 0.8780 0.420421
## Treatment 5.0455 0.009055 **
## InitialCommunity 132.8935 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment + InitialCommunity + (1 | Fragment) anova(LME_Ofav1.1) ranova(LME_Ofav1.1) # EMMs
Ofav.YII.emm1.1<-emmeans(LME_Ofav1.1, ~ InitialCommunity|Treatment)
# contrast(Ssid.YII.emm, "tukey")
# Tukey comparison (do not trust CI to compare EMMs)
plot(emmeans(LME_Ofav1.1, ~Treatment|InitialCommunity), comparisons = TRUE) +
coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_grid(~InitialCommunity) +
theme_bw()
Ofav.YII.emm1.1<-cld(Ofav.YII.emm1.1, by=NULL) # compact-letter display
Ofav.YII.emm1.1 #write.csv(Ofav.YII.emm1.1, "Outputs/Multicomp_OfavYII1.csv", row.names = F)- Model for “Days” as continous variable (days 1-76)
# 1. Ofav- Days
LME_Ofav<-lmer(YII ~ Treatment * Days * InitialCommunity +
(1|Genotype) + (1|Fragment) +(1|Replicate),
data=YII.Ofav.C, na.action=na.omit)
lmerTest::step (LME_Ofav) # Replicate is not significant## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 16 1243.0 -2454.1
## (1 | Replicate) 1 15 1243.0 -2456.1 0.000 1 1
## (1 | Genotype) 0 14 1233.7 -2439.4 18.642 1 1.577e-05 ***
## (1 | Fragment) 0 14 1233.5 -2438.9 19.149 1 1.209e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF
## Treatment:Days:InitialCommunity 0 0.0095086 0.0047543 2 570
## F value Pr(>F)
## Treatment:Days:InitialCommunity 5.128 0.006205 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment + Days + InitialCommunity + (1 | Genotype) +
## (1 | Fragment) + Treatment:Days + Treatment:InitialCommunity +
## Days:InitialCommunity + Treatment:Days:InitialCommunity ranova(LME_Ofav) LME_Ofav1<-lmer(YII ~ Treatment * Days * InitialCommunity +
(1 | Genotype/Fragment),
data=YII.Ofav.C, na.action=na.omit)
lmerTest::step (LME_Ofav1)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df
## <none> 15 1243.0 -2456.1
## (1 | Fragment:Genotype) 0 14 1233.5 -2438.9 19.149 1
## (1 | Genotype) 0 14 1233.7 -2439.4 18.642 1
## Pr(>Chisq)
## <none>
## (1 | Fragment:Genotype) 1.209e-05 ***
## (1 | Genotype) 1.577e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF
## Treatment:Days:InitialCommunity 0 0.0095086 0.0047543 2 570
## F value Pr(>F)
## Treatment:Days:InitialCommunity 5.128 0.006205 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment * Days * InitialCommunity + (1 | Genotype/Fragment)# 2. Predict values:
pred_Ofav1 <- predict(LME_Ofav1,re.form = NA)
#3. Bootstrap CI:
OF.boot1 <- bootMer(LME_Ofav1, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
std.err <- apply(OF.boot1$t, 2, sd)
CI.lo_1 <- pred_Ofav1 - std.err*1.96
CI.hi_1 <- pred_Ofav1 + std.err*1.96
#Plot
Model_of_1b_plot<- ggplot(
YII.Ofav.C, aes(x = Days, y = YII, colour = Treatment)) +
geom_line(aes(y = pred_Ofav1),size=2) +
#geom_point(aes(fill=factor(Treatment)),
# shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
size=2, alpha = 0.1, linetype = 0) +
#scale_color_manual(values=my_colours) +
#scale_fill_manual(values=my_colours) +
scale_y_continuous(name=expression(~italic("Fv / Fm")),
limits = c(0.35, 0.61),
breaks = seq(0.4, 0.6, by=0.1), expand = c(0,0))+
scale_x_continuous("Days in the experiment", limits = c(0, 78),
breaks = seq(0, 76, by=7), expand = c(0,0))+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position = position_dodge(1) )+
stat_summary(fun.y=mean, geom="line", position = position_dodge(1),
linetype=1, alpha=0.5) +
# stat_summary(fun.y=mean, geom="point", size =1,
# position=position_dodge(width=1), alpha=0.8) +
ggthe_bw
Model_of_1b_plot + facet_grid(~InitialCommunity)
- Model for “Days” as factor
# 1. Ofav- Days
LME_Ofav<-lmer(YII ~ Treatment * DaysF * InitialCommunity +
(1|Genotype) + (1|Fragment) +(1|Replicate),
data=YII.Ofav.C, na.action=na.omit)
lmerTest::step (LME_Ofav) # Replicate is not significant## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 58 1311.5 -2507.1
## (1 | Replicate) 1 57 1311.5 -2509.1 0.000 1 1
## (1 | Genotype) 0 56 1302.2 -2492.5 18.642 1 1.577e-05 ***
## (1 | Fragment) 0 56 1273.0 -2434.0 77.139 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF
## Treatment:DaysF:InitialCommunity 0 0.021644 0.0013528 16 528
## F value Pr(>F)
## Treatment:DaysF:InitialCommunity 2.7855 0.0002475 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment + DaysF + InitialCommunity + (1 | Genotype) +
## (1 | Fragment) + Treatment:DaysF + Treatment:InitialCommunity +
## DaysF:InitialCommunity + Treatment:DaysF:InitialCommunity ranova(LME_Ofav) LME_Ofav1.2<-lmer(YII ~ Treatment * DaysF * InitialCommunity +
(1|Fragment),
data=YII.Ofav.C, na.action=na.omit)
lmerTest::step (LME_Ofav1.2)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 56 1302.2 -2492.5
## (1 | Fragment) 0 55 1230.0 -2349.9 144.55 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF
## Treatment:DaysF:InitialCommunity 0 0.021644 0.0013528 16 528
## F value Pr(>F)
## Treatment:DaysF:InitialCommunity 2.7856 0.0002475 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment * DaysF * InitialCommunity + (1 | Fragment) # 2. Predict values:
pred_Ofav1.2 <- predict(LME_Ofav1.2,re.form = NA)
#3. Bootstrap CI:
OF.boot1.2 <- bootMer(LME_Ofav1.2, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
std.err <- apply(OF.boot1.2$t, 2, sd)
CI.lo_1 <- pred_Ofav1.2 - std.err*1.96
CI.hi_1 <- pred_Ofav1.2 + std.err*1.96
#Plot
Model_of_1.2_plot<- ggplot(
YII.Ofav.C, aes(x = Days, y = YII, colour = Treatment)) +
geom_line(aes(y = pred_Ofav1.2),size=2) +
#geom_point(aes(fill=factor(Treatment)),
# shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
size=2, alpha = 0.1, linetype = 0) +
#scale_color_manual(values=my_colours) +
#scale_fill_manual(values=my_colours) +
scale_y_continuous(name=expression(~italic("Fv / Fm")),
limits = c(0.35, 0.61),
breaks = seq(0.4, 0.6, by=0.1), expand = c(0,0))+
scale_x_continuous("Days in the experiment", limits = c(0, 78),
breaks = seq(0, 76, by=7), expand = c(0,0))+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position = position_dodge(1) )+
stat_summary(fun.y=mean, geom="point", size =1,
position=position_dodge(width=1), alpha=0.8) +
ggthe_bw
Model_of_1.2_plot + facet_grid(~InitialCommunity)
# EMMs
Ofav.YII.emmC<-emmeans(LME_Ofav1.2, ~Treatment * DaysF * InitialCommunity)
#contrast(Ofav.YII.emm, "tukey")
# Ofav.YII_groupsC<-cld(Ofav.YII.emmC, by=NULL) # compact-letter display
# Ofav.YII_groupsC<-Ofav.YII_groupsC[order(
# Ofav.YII_groupsC$Days,
# Ofav.YII_groupsC$Treatment,
# Ofav.YII_groupsC$InitialCommunity),]
# Ofav.YII_groupsC
# write.csv(Ofav.YII_groupsC, "Outputs/Multicomp_OfavYIIC.csv", row.names = F)Ssid
- Days pooled (8-76 Days)
LME_Ssid1.1<-lmer(YII ~ Treatment * InitialCommunity + (1|Fragment),
data=YII.Ssid.C.1)
step (LME_Ssid1.1) # Replicate is significant## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 2224.0 -4425.9
## (1 | Fragment) 0 10 2211.5 -4403.1 24.84 1 6.229e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF
## Treatment:InitialCommunity 0 0.039046 0.0097615 4 153.16
## F value Pr(>F)
## Treatment:InitialCommunity 5.9109 0.000189 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment * InitialCommunity + (1 | Fragment) anova(LME_Ssid1.1) ranova(LME_Ssid1.1) # EMMs
Ssid.YII.emm1.1<-emmeans(LME_Ssid1.1, ~ InitialCommunity|Treatment)
# contrast(Ssid.YII.emm, "tukey")
# Tukey comparison (do not trust CI to compare EMMs)
plot(emmeans(LME_Ssid1.1, ~Treatment|InitialCommunity), comparisons = TRUE) +
coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_grid(~InitialCommunity)+
theme_bw()
Ssid.YII.Groups1.1<-cld(Ssid.YII.emm1.1, by=NULL) # compact-letter display
Ssid.YII.Groups1.1<- Ssid.YII.Groups1.1[order(
Ssid.YII.Groups1.1$InitialCommunity,
Ssid.YII.Groups1.1$Treatment),]
Ssid.YII.Groups1.1 #write.csv( Ssid.YII.Groups1.1, "Outputs/Multicomp_SsidYII1.1.csv", row.names = F)- Model for “Days” as continous variable (days 1-76)
LME_Ssid<-lmer(YII ~ Treatment * Days * InitialCommunity +
(1|Fragment) +(1|Replicate),
data=subset(YII.Ssid.C))
step (LME_Ssid) # Replicate is not significant## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 21 2486.7 -4931.3
## (1 | Replicate) 1 20 2486.7 -4933.3 0.000 1 1
## (1 | Fragment) 0 19 2466.3 -4894.6 40.775 1 1.708e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF
## Treatment:Days:InitialCommunity 1 0.004115 0.0010287 4
## Treatment:Days 0 0.028806 0.0144029 2
## Treatment:InitialCommunity 0 0.033666 0.0084166 4
## Days:InitialCommunity 0 0.015350 0.0076750 2
## DenDF F value Pr(>F)
## Treatment:Days:InitialCommunity 1283.90 0.6706 0.6124001
## Treatment:Days 1288.52 9.3992 8.861e-05 ***
## Treatment:InitialCommunity 153.14 5.4926 0.0003694 ***
## Days:InitialCommunity 1287.90 5.0087 0.0068106 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment + Days + InitialCommunity + (1 | Fragment) +
## Treatment:Days + Treatment:InitialCommunity + Days:InitialCommunity anova(LME_Ssid) ranova(LME_Ssid) LME_Ssid1<-lmer(YII ~ Treatment * Days * InitialCommunity +
(1|Fragment),
data=YII.Ssid.C, na.action=na.omit)
# EMMs
Ssid.YII.emm1<-emmeans(LME_Ssid1, ~Treatment * Days * InitialCommunity)
# contrast(Ssid.YII.emm, "tukey")
# Effect plot options
emmip(LME_Ssid1, ~InitialCommunity|Treatment, CIs = TRUE) + theme_bw() # interaction plot of predictions
# 2. Predict values:
predSs1 <- predict(LME_Ssid1,re.form = NA)
# 3. Bootstrap CI:
boot1Ss <- bootMer(LME_Ssid1, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
std.err <- apply(boot1Ss$t, 2, sd)
CI.lo_1 <- predSs1 - std.err*1.96
CI.hi_1 <- predSs1 + std.err*1.96# 4. Plot
Model1Ss_plot<- ggplot(
YII.Ssid.C, aes(x = Days, y = YII, colour = Treatment)) +
geom_line(aes(y = predSs1),size=2) +
#geom_point(aes(fill=factor(Treatment)),
# shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.5) +
geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
size=2, alpha = 0.1, linetype = 0) +
#scale_color_manual(values=my_colours) +
#scale_fill_manual(values=my_colours) +
scale_y_continuous(name=expression(~italic("Fv / Fm")),
limits = c(0.3, 0.6),
breaks = seq(0.3, 0.6, by=0.1), expand = c(0,0))+
scale_x_continuous("Days in the experiment (C phase)", limits = c(0, 78),
breaks = seq(0, 76, by=7), expand = c(0,0))+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position = position_dodge(1) )+
stat_summary(fun.y=mean, geom="line", position = position_dodge(1),
linetype=1, alpha=0.5) +
#stat_summary(fun.y=mean, geom="point", size =1,
# position=position_dodge(width=1), alpha=0.8) +
ggthe_bw # + Fill.colour
Model1Ss_plot + facet_grid(~InitialCommunity)
- Model for “Days” as factor
# 1. Ssid- Days
LME_Ssid<-lmer(YII ~ Treatment * DaysF * InitialCommunity +
(1|Genotype) + (1|Fragment) +(1|Replicate),
data=YII.Ssid.C, na.action=na.omit)
lmerTest::step (LME_Ssid) # Replicate is not significant## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 85 3033.7 -5897.4
## (1 | Replicate) 1 84 3033.3 -5898.6 0.841 1 0.359
## (1 | Genotype) 0 83 2995.3 -5824.7 75.896 1 <2e-16 ***
## (1 | Fragment) 0 83 2972.7 -5779.4 121.188 1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF
## Treatment:DaysF:InitialCommunity 1 0.019444 0.0006076 32
## Treatment:DaysF 0 0.134201 0.0083875 16
## Treatment:InitialCommunity 0 0.020396 0.0050991 4
## DaysF:InitialCommunity 0 0.044362 0.0027727 16
## DenDF F value Pr(>F)
## Treatment:DaysF:InitialCommunity 1220.34 1.1718 0.2354
## Treatment:DaysF 1252.49 16.1060 < 2.2e-16 ***
## Treatment:InitialCommunity 147.69 9.7915 4.712e-07 ***
## DaysF:InitialCommunity 1252.33 5.3241 5.033e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment + DaysF + InitialCommunity + (1 | Genotype) +
## (1 | Fragment) + Treatment:DaysF + Treatment:InitialCommunity +
## DaysF:InitialCommunity ranova(LME_Ssid) LME_Ssid1.2<-lmer(YII ~ Treatment * DaysF * InitialCommunity +
(1|Fragment),
data=YII.Ssid.C, na.action=na.omit)
lmerTest::step (LME_Ssid1.2)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 83 2995.3 -5824.7
## (1 | Fragment) 0 82 2829.0 -5494.0 332.65 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF
## Treatment:DaysF:InitialCommunity 1 0.019389 0.0006059 32
## Treatment:DaysF 0 0.134162 0.0083851 16
## Treatment:InitialCommunity 0 0.011521 0.0028803 4
## DaysF:InitialCommunity 0 0.044336 0.0027710 16
## DenDF F value Pr(>F)
## Treatment:DaysF:InitialCommunity 1220.22 1.1685 0.2390861
## Treatment:DaysF 1252.32 16.1009 < 2.2e-16 ***
## Treatment:InitialCommunity 153.05 5.5307 0.0003475 ***
## DaysF:InitialCommunity 1252.23 5.3208 5.14e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment + DaysF + InitialCommunity + (1 | Fragment) +
## Treatment:DaysF + Treatment:InitialCommunity + DaysF:InitialCommunity# 2. Predict values:
pred_Ssid1.2 <- predict(LME_Ssid1.2,re.form = NA)
#3. Bootstrap CI:
Ss.boot1.2 <- bootMer(LME_Ssid1.2, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
std.err <- apply(Ss.boot1.2$t, 2, sd)
CI.lo_1 <- pred_Ssid1.2 - std.err*1.96
CI.hi_1 <- pred_Ssid1.2 + std.err*1.96# 4. Plot
Model_Ss_1.2_plot<- ggplot(
YII.Ssid.C, aes(x = Days, y = YII, colour = Treatment)) +
geom_line(aes(y = pred_Ssid1.2),size=2) +
#geom_point(aes(fill=factor(Treatment)),
# shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
size=2, alpha = 0.1, linetype = 0) +
#scale_color_manual(values=my_colours) +
#scale_fill_manual(values=my_colours) +
scale_y_continuous(name=expression(~italic("Fv / Fm")),
limits = c(0.35, 0.6),
breaks = seq(0.4, 0.6, by=0.1), expand = c(0,0))+
scale_x_continuous("Days in the experiment", limits = c(0, 78),
breaks = seq(0, 76, by=7), expand = c(0,0))+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position = position_dodge(1) )+
stat_summary(fun.y=mean, geom="point", size =1,
position=position_dodge(width=1), alpha=0.8) +
ggthe_bw
Model_Ss_1.2_plot + facet_grid(~InitialCommunity)
# 5. EMMs
Ssid.YII.emmC<-emmeans(LME_Ssid1.2, ~Treatment * DaysF * InitialCommunity)
#contrast(Ssid.YII.emm, "tukey")
Ssid.YII_groupsC<-cld(Ssid.YII.emmC, by=NULL) # compact-letter display
Ssid.YII_groupsC<-Ssid.YII_groupsC[order(
Ssid.YII_groupsC$DaysF,
Ssid.YII_groupsC$InitialCommunity,
Ssid.YII_groupsC$Treatment),]
Ssid.YII_groupsC #write.csv(Ssid.YII_groupsC, "Outputs/Multicomp_SsidYIIC.csv", row.names = F)H: Effect of pre-exposure to nutrient treatments during heat challenge
Ofav
- Model for “Days” as continous variable (days 84-110)
# 1. Ofav- Days
LME_Ofav.H<-lmer(YII ~ Treatment * Days * InitialCommunity +
(1|Genotype) + (1|Fragment) +(1|Replicate),
data=YII.Ofav.H, na.action=na.omit)
lmerTest::step (LME_Ofav.H) # Replicate is not significant## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 16 666.73 -1301.5
## (1 | Replicate) 1 15 666.73 -1303.5 0.0000 1 1.000000
## (1 | Genotype) 0 14 662.44 -1296.9 8.5781 1 0.003402 **
## (1 | Fragment) 0 14 663.16 -1298.3 7.1280 1 0.007589 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF
## Treatment:Days:InitialCommunity 1 0.0022744 0.0011372 2
## Treatment:InitialCommunity 2 0.0016208 0.0008104 2
## Treatment:Days 0 0.0189804 0.0094902 2
## Days:InitialCommunity 0 0.0088774 0.0088774 1
## DenDF F value Pr(>F)
## Treatment:Days:InitialCommunity 390.99 0.7378 0.478849
## Treatment:InitialCommunity 89.11 0.5261 0.592744
## Treatment:Days 390.47 6.1692 0.002302 **
## Days:InitialCommunity 393.25 5.7709 0.016757 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment + Days + InitialCommunity + (1 | Genotype) +
## (1 | Fragment) + Treatment:Days + Days:InitialCommunity ranova(LME_Ofav.H) LME_Ofav.H1<-lmer(YII ~ Treatment * Days +
(1 | Genotype/Fragment),
data=YII.Ofav.H, na.action=na.omit)
lmerTest::step (LME_Ofav.H1)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df
## <none> 9 687.23 -1356.5
## (1 | Fragment:Genotype) 0 8 683.59 -1351.2 7.2867 1
## (1 | Genotype) 0 8 674.23 -1332.5 26.0040 1
## Pr(>Chisq)
## <none>
## (1 | Fragment:Genotype) 0.006947 **
## (1 | Genotype) 3.407e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treatment:Days 0 0.015625 0.0078127 2 390 5.0081 0.00712
##
## Treatment:Days **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment * Days + (1 | Genotype/Fragment)# 2. Predict values:
pred_Ofav.H <- predict(LME_Ofav.H1,re.form = NA)
#3. Bootstrap CI:
OF.bootH <- bootMer(LME_Ofav.H1, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
std.err <- apply(OF.bootH$t, 2, sd)
CI.lo_H <- pred_Ofav.H - std.err*1.96
CI.hi_H <- pred_Ofav.H + std.err*1.96 # 4 .Plot
Model_of_H_plot<- ggplot(
YII.Ofav.H, aes(x = Days, y = YII, colour = Treatment)) +
geom_line(aes(y = pred_Ofav.H),size=2) +
#geom_point(aes(fill=factor(Treatment)),
# shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
geom_ribbon(aes(ymin = CI.lo_H, ymax = CI.hi_H),
size=2, alpha = 0.1, linetype = 0) +
#scale_color_manual(values=my_colours) +
#scale_fill_manual(values=my_colours) +
scale_y_continuous(name=expression(~italic("Fv / Fm")),
limits = c(0.25, 0.58),
breaks = seq(0.3, 0.5, by=0.1), expand = c(0,0))+
scale_x_continuous("Days in the experiment (H phase)", limits = c(80, 112),
breaks = seq(80, 110, by=7), expand = c(0,0))+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position = position_dodge(1) )+
stat_summary(fun.y=mean, geom="line", position = position_dodge(1),
linetype=1, alpha=1) +
stat_summary(fun.y=mean, geom="point", size =1,
position=position_dodge(width=1), alpha=0.8) +
ggthe_bw
Model_of_H_plot
Model_of_H_plot+ facet_grid(InitialCommunity~.)
- Model for “Days” as factor (daysF 84-110)
# 1. Ofav- Days
LME_Ofav.HF<-lmer(YII ~ Treatment * DaysF * Community +
(1|Genotype) + (1|Fragment) +(1|Replicate),
data=YII.Ofav.H, na.action=na.omit)
lmerTest::step (LME_Ofav.HF) # Replicate is not significant## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 49 740.58 -1383.2
## (1 | Replicate) 1 48 740.58 -1385.2 0.000 1 1.000000
## (1 | Genotype) 0 47 735.75 -1377.5 9.650 1 0.001894 **
## (1 | Fragment) 0 47 715.64 -1337.3 49.865 1 1.647e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF
## Treatment:DaysF:Community 1 0.001743 0.0002179 8 330.04
## Treatment:Community 2 0.000134 0.0000672 2 110.06
## Treatment:DaysF 0 0.033700 0.0021062 16 330.20
## DaysF:Community 0 0.047109 0.0067299 7 336.98
## F value Pr(>F)
## Treatment:DaysF:Community 0.3417 0.9492
## Treatment:Community 0.1071 0.8986
## Treatment:DaysF 3.3619 1.663e-05 ***
## DaysF:Community 10.7421 3.153e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment + DaysF + Community + (1 | Genotype) + (1 | Fragment) +
## Treatment:DaysF + DaysF:Community ranova(LME_Ofav.HF) LME_Ofav.HF.1<-lmer(YII ~ Treatment + DaysF + Community +
(1 | Genotype) + (1 | Fragment) +
Treatment:DaysF + DaysF:Community,
data=YII.Ofav.H, na.action=na.omit)
lmerTest::step (LME_Ofav.HF.1)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 38 767.61 -1459.2
## (1 | Genotype) 0 37 762.53 -1451.1 10.158 1 0.001437 **
## (1 | Fragment) 0 37 742.39 -1410.8 50.450 1 1.222e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## Treatment:DaysF 0 0.033700 0.0021062 16 330.20 3.3619
## DaysF:Community 0 0.047109 0.0067299 7 336.98 10.7421
## Pr(>F)
## Treatment:DaysF 1.663e-05 ***
## DaysF:Community 3.153e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment + DaysF + Community + (1 | Genotype) + (1 | Fragment) +
## Treatment:DaysF + DaysF:Community# 2. Predict values:
pred_OfavHF <- predict(LME_Ofav.HF.1,re.form = NA)
# 3. Bootstrap CI:
OF.bootHF <- bootMer(LME_Ofav.HF.1, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
std.err <- apply(OF.bootHF$t, 2, sd)
CI.lo_HF <- pred_OfavHF - std.err*1.96
CI.hi_HF <- pred_OfavHF + std.err*1.96# 4. Plot
Model_of_HF_plot<- ggplot(
YII.Ofav.H, aes(x = Days, y = YII, colour = Treatment)) +
geom_line(aes(y = pred_OfavHF),size=2) +
#geom_point(aes(fill=factor(Treatment)),
# shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
geom_ribbon(aes(ymin = CI.lo_HF, ymax = CI.hi_HF),
size=2, alpha = 0.1, linetype = 0) +
#scale_color_manual(values=my_colours) +
#scale_fill_manual(values=my_colours) +
scale_y_continuous(name=expression(~italic("Fv / Fm")),
limits = c(0.25, 0.55),
breaks = seq(0.3, 0.5, by=0.1), expand = c(0,0))+
scale_x_continuous("Days in the experiment (H phase)", limits = c(80, 112),
breaks = seq(80, 110, by=7), expand = c(0,0))+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position = position_dodge(1) )+
stat_summary(fun.y=mean, geom="point", size =1,
position=position_dodge(width=1), alpha=0.8) +
ggthe_bw
Model_of_HF_plot
Model_of_HF_plot + facet_grid(Community~.)
# 5. EMMs
Ofav.YII.emmH<-emmeans(LME_Ofav.HF.1, ~ Treatment:DaysF + DaysF:Community)
#contrast(Ofav.YII.emmH, "tukey")
Ofav.YII_groupsH<-cld(Ofav.YII.emmH, by=NULL) # compact-letter display
Ofav.YII_groupsH<-Ofav.YII_groupsH[order(
Ofav.YII_groupsH$Days,
Ofav.YII_groupsH$Treatment,
Ofav.YII_groupsH$Community),]
Ofav.YII_groupsH #write.csv(Ofav.YII_groupsH, "Outputs/Multicomp_OfavYII_H.csv", row.names = F)Ssid
- Model for “Days” as continous variable (days 84-110)
# 1. Ssid- Days
LME_Ssid.H<-lmer(YII ~ Treatment * Days * InitialCommunity +
(1|Genotype) + (1|Fragment) +(1|Replicate),
data=YII.Ssid.H, na.action=na.omit)
lmerTest::step (LME_Ssid.H) # Replicate is not significant## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 22 1692.9 -3341.9
## (1 | Replicate) 1 21 1692.9 -3343.9 0.0000 1 1.000000
## (1 | Genotype) 0 20 1690.8 -3341.5 4.3147 1 0.037785 *
## (1 | Fragment) 0 20 1688.9 -3337.9 7.9765 1 0.004739 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF
## Treatment:Days:InitialCommunity 0 0.043395 0.010849 4 1029
## F value Pr(>F)
## Treatment:Days:InitialCommunity 4.6856 0.0009452 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment + Days + InitialCommunity + (1 | Genotype) +
## (1 | Fragment) + Treatment:Days + Treatment:InitialCommunity +
## Days:InitialCommunity + Treatment:Days:InitialCommunity ranova(LME_Ssid.H) LME_Ssid.H1<-lmer(YII ~ Treatment * Days * InitialCommunity +
(1 | Genotype/Fragment),
data=YII.Ssid.H, na.action=na.omit)
lmerTest::step (LME_Ssid.H1)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df
## <none> 21 1692.9 -3343.9
## (1 | Fragment:Genotype) 0 20 1688.9 -3337.9 7.9765 1
## (1 | Genotype) 0 20 1690.8 -3341.5 4.3147 1
## Pr(>Chisq)
## <none>
## (1 | Fragment:Genotype) 0.004739 **
## (1 | Genotype) 0.037785 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF
## Treatment:Days:InitialCommunity 0 0.043395 0.010849 4 1029
## F value Pr(>F)
## Treatment:Days:InitialCommunity 4.6856 0.0009452 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment * Days * InitialCommunity + (1 | Genotype/Fragment) anova(LME_Ssid.H1)# 2. Predict values:
pred_Ssid.H <- predict(LME_Ssid.H1,re.form = NA)
#3. Bootstrap CI:
Ss.bootH <- bootMer(LME_Ssid.H1, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
std.errSH <- apply(Ss.bootH$t, 2, sd)
CI.lo_SH <- pred_Ssid.H - std.errSH*1.96
CI.hi_SH <- pred_Ssid.H + std.errSH*1.96 # 4 .Plot
Model_Ss_H_plot<- ggplot(
YII.Ssid.H, aes(x = Days, y = YII, colour = Treatment)) +
geom_line(aes(y = pred_Ssid.H),size=2) +
#geom_point(aes(fill=factor(Treatment)),
# shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
geom_ribbon(aes(ymin = CI.lo_SH, ymax = CI.hi_SH),
size=2, alpha = 0.1, linetype = 0) +
#scale_color_manual(values=my_colours) +
#scale_fill_manual(values=my_colours) +
scale_y_continuous(name=expression(~italic("Fv / Fm")),
limits = c(0.1, 0.58),
breaks = seq(0.1, 0.5, by=0.1), expand = c(0,0))+
scale_x_continuous("Days in the experiment (H phase)", limits = c(80, 112),
breaks = seq(80, 110, by=7), expand = c(0,0))+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position = position_dodge(1) )+
stat_summary(fun.y=mean, geom="line", position = position_dodge(1),
linetype=1, alpha=0.5) +
#stat_summary(fun.y=mean, geom="point", size =1,
# position=position_dodge(width=1), alpha=0.8) +
ggthe_bw
Model_Ss_H_plot+ facet_grid(~InitialCommunity)
- Model for “Days” as factor (daysF 84-110)
# 1. Ssid- Days
LME_Ssid.HF<-lmer(YII ~ Treatment * DaysF * InitialCommunity +
(1|Genotype) + (1|Fragment) +(1|Replicate),
data=YII.Ssid.H, na.action=na.omit)
lmerTest::step (LME_Ssid.HF) # Replicate is not significant## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 85 1885.7 -3601.5
## (1 | Replicate) 1 84 1885.7 -3603.5 -0.010 1 1.0000
## (1 | Genotype) 0 83 1883.2 -3600.3 5.153 1 0.0232 *
## (1 | Fragment) 0 83 1846.5 -3526.9 78.534 1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF
## Treatment:DaysF:InitialCommunity 0 0.081998 0.0025624 32
## DenDF F value Pr(>F)
## Treatment:DaysF:InitialCommunity 929.23 2.275 8.27e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment + DaysF + InitialCommunity + (1 | Genotype) +
## (1 | Fragment) + Treatment:DaysF + Treatment:InitialCommunity +
## DaysF:InitialCommunity + Treatment:DaysF:InitialCommunity ranova(LME_Ssid.HF) LME_Ssid.HF.1<-lmer(YII ~ Treatment * DaysF * InitialCommunity + (1 | Fragment),
data=YII.Ssid.H, na.action=na.omit)
lmerTest::step (LME_Ssid.HF.1)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 83 1883.2 -3600.3
## (1 | Fragment) 0 82 1835.2 -3506.5 95.807 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF
## Treatment:DaysF:InitialCommunity 0 0.081916 0.0025599 32
## DenDF F value Pr(>F)
## Treatment:DaysF:InitialCommunity 928.77 2.2693 8.706e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## YII ~ Treatment * DaysF * InitialCommunity + (1 | Fragment) anova(LME_Ssid.HF.1)# 2. Predict values:
pred_SsidHF <- predict(LME_Ssid.HF.1,re.form = NA)
# 3. Bootstrap CI:
Ss.bootHF <- bootMer(LME_Ssid.HF.1, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
std.err <- apply(Ss.bootHF$t, 2, sd)
CI.lo_HF.S <- pred_SsidHF - std.err*1.96
CI.hi_HF.S <- pred_SsidHF + std.err*1.96# 4. Plot
Model_Ss_HF_plot<- ggplot(
YII.Ssid.H, aes(x = Days, y = YII, colour = Treatment)) +
geom_line(aes(y = pred_SsidHF),size=2) +
#geom_point(aes(fill=factor(Treatment)),
# shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
geom_ribbon(aes(ymin = CI.lo_HF.S, ymax = CI.hi_HF.S),
size=2, alpha = 0.1, linetype = 0) +
#scale_color_manual(values=my_colours) +
#scale_fill_manual(values=my_colours) +
scale_y_continuous(name=expression(~italic("Fv / Fm")),
limits = c(0.1, 0.55),
breaks = seq(0.1, 0.5, by=0.1), expand = c(0,0))+
scale_x_continuous("Days in the experiment (H phase)", limits = c(80, 112),
breaks = seq(80, 110, by=7), expand = c(0,0))+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position = position_dodge(1) )+
#stat_summary(fun.y=mean, geom="point", size =1,
# position=position_dodge(width=1), alpha=0.8) +
ggthe_bw
Model_Ss_HF_plot + facet_grid(~InitialCommunity)
# 5. EMMs
Ssid.YII.emmH<-emmeans(LME_Ssid.HF.1, ~ Treatment*DaysF*InitialCommunity)
#contrast(Ssid.YII.emmH, "tukey")
Ssid.YII_groupsH<-cld(Ssid.YII.emmH, by=NULL) # compact-letter display
Ssid.YII_groupsH<-Ssid.YII_groupsH[order(
Ssid.YII_groupsH$Days,
Ssid.YII_groupsH$Treatment,
Ssid.YII_groupsH$InitialCommunity),]
Ssid.YII_groupsH #write.csv(Ssid.YII_groupsH, "Outputs/Multicomp_SsidYIIH.csv", row.names = F)Packages used
# Creates bibliography
#knitr::write_bib(c(.packages()), "packages.bib")Arnold, Jeffrey B. 2019. Ggthemes: Extra Themes, Scales and Geoms for ’Ggplot2’. https://CRAN.R-project.org/package=ggthemes.
Bates, Douglas, and Martin Maechler. 2019. Matrix: Sparse and Dense Matrix Classes and Methods. https://CRAN.R-project.org/package=Matrix.
Bates, Douglas, Martin Maechler, Ben Bolker, and Steven Walker. 2019. Lme4: Linear Mixed-Effects Models Using ’Eigen’ and S4. https://CRAN.R-project.org/package=lme4.
Fox, John, Sanford Weisberg, and Brad Price. 2018. CarData: Companion to Applied Regression Data Sets. https://CRAN.R-project.org/package=carData.
Fox, John, Sanford Weisberg, Brad Price, Michael Friendly, and Jangman Hong. 2019. Effects: Effect Displays for Linear, Generalized Linear, and Other Models. https://CRAN.R-project.org/package=effects.
Genz, Alan, Frank Bretz, Tetsuhisa Miwa, Xuefei Mi, and Torsten Hothorn. 2019. Mvtnorm: Multivariate Normal and T Distributions. https://CRAN.R-project.org/package=mvtnorm.
Gohel, David, Hadley Wickham, Lionel Henry, and Jeroen Ooms. 2019. Gdtools: Utilities for Graphical Rendering. https://CRAN.R-project.org/package=gdtools.
Henry, Lionel, and Hadley Wickham. 2019. Purrr: Functional Programming Tools. https://CRAN.R-project.org/package=purrr.
Hothorn, Torsten. 2019. TH.data: TH’s Data Archive. https://CRAN.R-project.org/package=TH.data.
Hothorn, Torsten, Frank Bretz, and Peter Westfall. 2019. Multcomp: Simultaneous Inference in General Parametric Models. https://CRAN.R-project.org/package=multcomp.
Kuznetsova, Alexandra, Per Bruun Brockhoff, and Rune Haubo Bojesen Christensen. 2019. LmerTest: Tests in Linear Mixed Effects Models. https://CRAN.R-project.org/package=lmerTest.
Lenth, Russell. 2019. Emmeans: Estimated Marginal Means, Aka Least-Squares Means. https://CRAN.R-project.org/package=emmeans.
Müller, Kirill, and Hadley Wickham. 2019. Tibble: Simple Data Frames. https://CRAN.R-project.org/package=tibble.
R Core Team. 2020. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Ripley, Brian. 2019. MASS: Support Functions and Datasets for Venables and Ripley’s Mass. https://CRAN.R-project.org/package=MASS.
Therneau, Terry M. 2019. Survival: Survival Analysis. https://CRAN.R-project.org/package=survival.
Wickham, Hadley. 2017. Tidyverse: Easily Install and Load the ’Tidyverse’. https://CRAN.R-project.org/package=tidyverse.
———. 2019a. Forcats: Tools for Working with Categorical Variables (Factors). https://CRAN.R-project.org/package=forcats.
———. 2019b. Stringr: Simple, Consistent Wrappers for Common String Operations. https://CRAN.R-project.org/package=stringr.
Wickham, Hadley, and Lionel Henry. 2020. Tidyr: Tidy Messy Data. https://CRAN.R-project.org/package=tidyr.
Wickham, Hadley, Winston Chang, Lionel Henry, Thomas Lin Pedersen, Kohske Takahashi, Claus Wilke, Kara Woo, and Hiroaki Yutani. 2019. Ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics. https://CRAN.R-project.org/package=ggplot2.
Wickham, Hadley, Romain François, Lionel Henry, and Kirill Müller. 2019. Dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.
Wickham, Hadley, Jim Hester, and Romain Francois. 2018. Readr: Read Rectangular Text Data. https://CRAN.R-project.org/package=readr.